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Water  management  plays  an  important  role  in  the  durability  and  efficiency  of  a  proton  exchange  mem¬ 
brane  fuel  cell  (PEMFC).  In  this  study,  a  single  cell  is  modeled  as  a  lumped  model  consisting  of  15 
interconnected  segments,  which  are  linked  according  to  the  flow  field  patterns  of  the  anode  and  cathode 
but  they  are  treated  as  lumped  elements  individually.  Parameters  of  this  model  were  calibrated  based  on 
neutron  radiography  experimental  results  obtained  at  the  NIST  Center  for  Neutron  Research  (NCNR).  The 
model  can  be  used  to  predict  distributions  of  current  density,  water  content  in  the  membrane,  relative 
humidity  (RH)  in  the  flow  channels,  and  water  accumulation  in  the  gas  diffusion  layer  (GDL). 

©  2008  Elsevier  B.V.  All  rights  reserved. 


1.  Introduction 

Many  factors  influence  the  performance  of  proton  exchange 
membrane  fuel  cells  (PEMFCs),  including  membrane  material  and 
thickness,  platinum  loading,  flow  field  designs,  temperature,  reac¬ 
tant  partial  pressure,  etc.  In  addition,  a  critical  issue  that  has 
severely  limited  the  application  of  PEMFCs  is  its  poor  reliability 
under  cyclic  temperature  and  humidity  operations.  An  important 
factor  that  influences  both  the  nominal  performance  and  life  under 
transient  loading  is  the  water  accumulation  and  distribution  in  a 
PEMFC.  When  the  fuel  cell  has  too  little  or  too  much  water,  both 
performance  and  reliability  suffers.  Water  accumulation  also  influ¬ 
ences  the  warm-up  and  shutdown  procedure  for  PEMFCs  that  need 
to  work  in  below-freezing  temperature.  It  is  fair  to  say  that  in 
addition  to  the  cost  issue,  water  management  is  one  of  the  most 
important  remaining  issues  for  the  adoption  of  PEMFCs.  Many 
models  have  been  developed  over  the  past  several  years  for  the 
water/humidity  behavior  inside  a  PEMFC. 

Bernardi  and  Verbrugge  [1,2]  developed  one  of  the  early 
steady-state,  one-dimensional  mathematical  models.  Their  model 
describes  reactant  transport  in  the  GDLs  and  water  transport  in  a 
PEMFC.  The  membrane  was  assumed  to  be  fully  hydrated,  which  is 
different  from  practical  working  conditions  of  PEMFCs,  especially 
for  the  anode  side.  Springer  et  al.  [3]  empirically  related  membrane 
conductivity  to  water  content  of  the  Nafion  membrane.  Many  of 
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the  subsequently  developed  models  used  this  empirical  relation¬ 
ship  to  determine  the  conductivity  of  the  Nafion  membrane,  even 
for  membranes  with  different  thickness.  Fuller  and  Newman  [4] 
developed  a  two-dimensional  model  to  discuss  water  management, 
thermal  management  and  fuel  utilization  in  a  PEMFC.  Gurau  et  al. 
[5]  developed  a  two-dimensional  model,  which  discussed  reactant 
concentrations  in  the  through-MEA  direction  and  along  the  flow 
channel  direction.  Um  et  al.  [6]  and  Wang  et  al.  [7]  developed  mod¬ 
els  based  on  computational  fluid  dynamics  (CFD)  and  solved  the 
equations  numerically.  Um  et  al.  [8]  also  extended  their  work  to  a 
three-dimensional  model  to  study  the  performance  of  an  interdig- 
itated  flow  field  design.  Their  results  show  that  forced  convection 
of  gases  through  GDL  helps  to  improve  performance  at  high  cur¬ 
rent  densities.  None  of  the  models  discussed  above  considered  the 
effect  of  water  accumulation  on  cell  performance. 

Starting  from  around  the  turn  of  the  century,  models  that 
include  the  water/humidity  behavior  start  to  appear  in  the  litera¬ 
ture.  Baschuk  et  al.  [9]  developed  a  model  with  the  effect  of  variable 
degree  of  water  flooding  in  the  cathode  catalyst  layer  and  cathode 
GDL  on  cell  performance.  Wang  et  al.  [10]  of  the  Pennsylvania  State 
University  developed  a  model  which  handles  the  situation  where 
two-phase  flow  exists  in  the  cathode.  Pasaogullari  and  Wang  [11] 
applied  the  two-phase  flow  model  in  the  cathode  GDL  and  investi¬ 
gated  the  effect  of  liquid  saturation  on  cell  performance.  Wang  et  al. 
[10]  and  Pasaogullari  and  Wang  [11  ]  related  capillary  pressure  with 
the  Leverett’s  function.  The  two-phase  flow  model  successfully 
described  water  vapor  distribution  and  liquid  water  accumulation 
in  the  GDL  and  in  the  flow  channel.  However,  most  of  these  models 
focus  on  water  accumulation  in  the  GDL  under  the  channel  and  the 


0378-7753 /$  -  see  front  matter  ©  2008  Elsevier  B.V.  All  rights  reserved. 
doi:10.1016/j.jpowsour.2008.07.018 


1180 


Y.-S.  Chen,  H.  Peng  /  Journal  of  Power  Sources  185  (2008)  1179-1192 


Nomenclature 

a  water  vapor  activity 

A  area(m2) 

c  concentration  (mol  m-3 ) 

d  hydraulic  diameter  (m) 

Di_j  diffusivity  of  gas  pair  i-j  in  a  mixture  (m2  s-1 ) 

F  Faraday’s  constant  (96485  °C  equiv.-1 ) 

H  channel  depth  (m) 

i  current  density  (A  m-2 ) 

I  current  (A) 

^osmotic  electro-osmotic  drag  coefficient 
Kdi ff  back  diffusion  coefficient  (mol  s  m-2 ) 

I<conv  coefficient  of  convective  mass  transfer  (mol  s  nrr2 ) 

L  channel  length  (m) 

M  equivalent  weight  of  a  dry  membrane  (kg  mol-1 ) 

N  molar  flow  rate  (mol  s_1) 

P  pressure  (Pa) 

Q  gas  volume  flow  rate  (m3  s-1 ) 

R  universal  gas  constant  (8.314J  mol-1  K-1 ) 

Rj  resistance  of  component  j  (£2) 

Sh  Sherwood  number 

t  thickness  (m) 

T  temperature  (I<) 

V  voltage  (V) 

W  channel  or  rib  width  (m) 

Xj  molar  fraction  of  species  j 

yo2  percentage  of  oxygen  in  the  air 

Z  channel  number  in  a  segment 

Greek  letters 

ac  transfer  coefficient 

s  porosity  of  gas  diffusion  layer 

cp  relative  humidity 

X  water  content 

fi  dynamic  viscosity  (kg  nrr1  s-1 ) 

p  density  (kg  nr3) 

a  electrical  conductivity  (£2_1  nr1 ) 

f  stoichiometry  of  gas 

Subscripts 
an  anode 

act  activation 

avg  average 

c  critical 

ca  cathode 

cap  capillary 

c/g  channel  and  gas  diffusion  layer  interface 

ch  channel 

contact  contact 

cone  concentration 

gdl  gas  diffusion  layer 

g/m  gas  diffusion  layer  and  membrane  interface 

H2  hydrogen 

H20  water 

in  inlet 

limit  limit 

N2  nitrogen 

ohm  ohmic 

out  outlet 

02  oxygen 

p  pore 


pern  proton  exchange  membrane 
plate  plate 

sat  saturation 

seg  segment 

v  water  vapor 

w  liquid  water 


developed  formula  do  not  apply  readily  to  the  subspace  under  the 
rib. 

Natarajan  and  Nguyen  [12]  from  the  University  of  Kansas  pro¬ 
posed  a  model  that  included  the  effect  of  water  accumulation  in  the 
GDL  under  the  rib  and  under  the  channel  on  cell  performance.  In 
their  model,  instead  of  using  the  Leverett’s  function,  they  suggested 
another  empirical  equation  to  describe  the  capillary  pressure.  Later, 
the  same  group  [13,14]  further  simplified  capillary  pressure  gradi¬ 
ent  in  their  models  as  a  constant.  Their  results  showed  significant 
difference  with  studies  using  the  Leverett’s  function.  Thus,  exper¬ 
imental  data  that  clearly  describe  water  accumulation  in  the  GDL, 
for  both  along  the  flow  direction  and  across  the  GDL  direction 
is  needed.  Recently,  several  studies  used  neutron  radiography  to 
detect  liquid  water  distribution  in  PEMFCs  [15-17].  Neutron  image 
results  showed  that  significant  amount  of  liquid  water  could  accu¬ 
mulate  in  the  GDL  under  the  ribs.  Therefore,  it  is  important  to 
develop  a  model  that  predicts  water  distribution  in  GDL  under  the 
rib  as  well  as  under  the  flow  channel. 

Because  the  reactant  concentration  varies  along  the  flow  chan¬ 
nels,  it  causes  variations  in  current  density,  water  content,  and 
temperature  [18-20].  Therefore,  the  water  generation  and  distri¬ 
bution  in  a  PEMFC  are  not  uniform.  In  addition,  different  anode  and 
cathode  flow  field  patterns  were  designed  for  different  applications 
or  working  conditions  [21-23].  Many  CFD  models  have  difficulties 
simulating  PEMFCs  with  complex  flow  fields  due  to  requirement  of 
heavy  numerical  computation  load.  Currently  published  CFD  mod¬ 
els  simulate  the  reaction  either  in  a  straight  flow  channel  or  in  a 
simple  flow  field.  Lumped  models  [24-26]  commonly  assume  a 
uniform  reaction  within  fuel  cells  and  do  not  consider  the  spa¬ 
tial  distribution  of  reactants.  Therefore,  pure  CFD  models  or  pure 
lumped  models  may  not  be  the  best  modeling  choice. 

In  this  study,  a  steady-state,  segmented  mathematical  model 
was  developed  to  describe  distributions  of  liquid  water  accumula¬ 
tion,  current  density,  and  relative  humidity  (RFI)  in  the  flow  channel 
of  a  PEMFC.  This  model  was  calibrated  by  using  neutron  radiogra¬ 
phy  experiments  to  quantify  liquid  water  in  a  PEMFC  with  the  same 
flow  field  pattern.  Water  transport  in  the  MEA  and  the  influence  of 
RH  of  cathode  inlet  is  also  discussed  in  this  study. 

2.  Mathematical  model 

To  capture  distributed  characteristics  of  a  PEMFC,  the  active 
area  is  divided  into  15  segments  that  are  connected  according  to 
flow  fields,  as  shown  in  Fig.  1.  Each  segment  is  viewed  as  a  small 
lumped  model,  i.e.  reactant/membrane  properties  and  reaction  in 
each  segment  are  assumed  to  be  uniform.  The  segments  are  con¬ 
nected  together  based  on  the  flow  direction  of  the  reactants.  Since 
each  is  regarded  as  a  lumped  model,  it  cannot  account  for  the  rib 
effects  on  gas  transport.  Flowever,  we  will  introduce  semi-empirical 
correlation  to  describe  the  rib  effects  on  water  accumulation. 

The  inputs  of  a  segment  are  the  outputs  of  the  preceding  seg¬ 
ments.  For  the  overall  cell,  input  variables  are  stoichiometry  value, 
RH,  and  temperatures  of  the  inflow  gas  and  cell  temperature. 
According  to  the  experiment  of  Wang  et  al.  [27],  temperature  differ¬ 
ence  between  upper  stream  and  down  stream  is  less  than  2  °C  when 
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Fig.  1.  Schematic  of  a  single  cell  modeled  as  several  small  segments. 


cell  current  density  is  0.74  A  cm-2.  The  study  of  Chen  and  Peng  [28] 
also  shows  uniformly  distributed  temperature  in  a  single  cell;  thus, 
all  segment  temperatures  are  assumed  the  same  and  constant  as 
operating  temperature.  The  assumptions  are  summarized  below: 


1.  The  model  describes  steady-state  conditions. 

2.  The  ideal  gas  law  was  employed  for  gas  mixture. 

3.  Temperature  throughout  the  single  cell  is  uniformly  distributed. 

4.  Chemical  reaction  throughout  the  segment  is  uniform. 

5.  Rib  effects  on  gas  transport  and  gas  transport  in  the  flow  chan¬ 
nel  direction  are  neglected.  Only  gas  transport  through  the  MEA 
direction  is  considered  in  each  segment. 


Based  on  the  desired  cell  current  and  operating  conditions,  the 
molar  flow  rates  of  inflow  hydrogen  and  oxygen  for  a  single  cell  are 
evaluated  as 


/ u  _  ^cell  v 

iVan,H2,in  —  2 p  ^ai 


M  —  icell  v. 

ivca,02,in  -  ^  S< 


(1) 

(2) 


where  N  is  the  molar  flow  rate  in  mols-1,  £an  and  £Ca  are  stoi¬ 
chiometry  of  anode  and  cathode,  respectively,  and  F  is  the  Faraday 
constant. 

If  air  is  used  as  the  cathode  reactant,  the  nitrogen  molar  flow 
rate  is  calculated  from 


based  on  the  flow  fields  of  anode  and  cathode  to  the  subsequent 
segment. 

Species  flow  in  each  segment  is  shown  in  Fig.  2.  Each  segment 
itself  consists  of  six  interacting  sub-models:  cathode  flow  channel, 
anode  flow  channel,  cathode  GDL,  anode  GDL,  membrane  hydra¬ 
tion,  and  segment  voltage.  Theses  models  will  be  described  in  the 
following  sections. 


2.1.  Anode/cathode  channel  model 

The  channel  model  describes  the  reactant  behavior  inside  the 
anode  and  cathode  of  a  segment.  The  model  uses  the  molar  con¬ 
servation  principle  and  fluid  dynamic  properties  to  calculate  the 
outflow  properties  and  pressure  drop  along  the  flow  channels.  The 
pressure  drop  of  the  gas  mixture  in  the  flow  channels  was  fre¬ 
quently  ignored  in  earlier  models;  however,  in  practice  it  is  one 
of  the  key  parameters  in  designing  a  fuel  cell  and  is  related  to  the 
selection  of  the  air  pump  and  the  calculation  of  efficiency  of  a  fuel 
cell  system. 

The  segment  current,  Jse g,  is  an  input  based  on  which  the  seg¬ 
ment  model  can  be  simulated.  The  amount  of  consumed  reactants 
in  the  segment  can  be  determined  by 

Nan,H2,  react  =  "ypr  (9) 


^ca,N2,in 


—  ^ca,02,in 


i  -yo2 

yo2 


(3) 


where  yo2  is  the  percentage  of  oxygen  in  the  air.  For  an  operating 
condition  at  selected  inlet  relative  humidity  cp  and  pressure  P,  the 
molar  fraction  of  inlet  vapor  can  be  calculated  from 


*an,v,in 


*ca,v,in 


^an,in^v,sat 

P  • 

1  an, in 
^ca,in^v,sat 

P  ■ 

1  ca,in 


(4) 

(5) 


where  Pv,Sat  is  the  saturated  vapor  pressure,  which  is  a  function 
of  working  temperature.  The  value  of  Pv,Sat  is  calculated  from  the 
following  equation,  given  in  Ref.  [3]: 

Pv,  sat  =  1.013  X  105 


xlQ[-2.1794+0.02953Tseg-9.1837xl0-5Ts2eg+1.4454xl0-7rs3eg] 


(6) 


The  inlet  water  molar  flow  rate  can  then  be  calculated  from 


Nai 


1  -xai 


-N; 


an,H2,in 


at  —  /Vca,v,in 

^ca,w,in  — 

1  *ca,v,i 


'(^ca,02,in  +  i^ca,N2,in) 


(7) 

(8) 


Eqs.  (l)-(8)  describe  required  amount  of  inflow  species.  After  fed 
into  the  first  segment  of  a  fuel  cell,  gases  flow  through  each  segment 


Nca,02,  react 


heg 

4  F 


(10) 


Anode 

Anode 

PEM 

Cathode 

Cathode 

Channel 

GDL 

GDL 

Channel 

No2,out 
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The  amount  of  water  generation  at  the  cathode  catalyst  layer 
can  be  expressed  as 


N  -  /seg 

^ca,w,gen  — 

(u) 

Then  the  molar  outflow  rate  can  be  calculated  from 

Nan, H2, out  —  Nan?H2,in  —  Nan?H2, react 

(12) 

NCa,02,out  —  Nca  o2,in  —  Nca?o2, react 

(13) 

Nan, w, out  —  Nanwin  —  Nan?wgdi 

(14) 

NCa,w,out  =  Nca?wjn  —  Nca?w?gdi 

(15) 

The  water  molar  flow  rates  Nan  wgdi 

and  Nca  wgdl  through  GDL  in 

Eqs.  (14)  and  (15)  are  calculated  by  the  gas  diffusion  layer  (GDL) 
model. 

The  pressure  drop  due  to  friction  for  a  continuous,  straight  chan¬ 
nel  with  length  L  can  be  calculated  from  the  following  equation 
[29]: 

ap  =  32  [  t5Y}Q{y)d  (16) 

Jo  ZAchdl 

where  the  reactant  is  consumed  along  the  flow  channel  direction  y. 
As  a  first-order  approximation,  we  assume  that  dynamic  viscosity 
/ji(y)  and  flow  rate  Q(y)  of  gas  mixtures  vary  linearly  along  the  flow 
channel: 

My)  —  M  in  —  J^iP'in  —  P'  out)  (17) 

Q(y)  =  Qin  -  l (Qin  -  Qout)  (18) 

where  /xin  and  /xout  are  the  dynamic  viscosities  of  the  mixture  at  the 
inlet  and  outlet,  respectively,  and  can  be  calculated  by  the  mixture 
properties 


n  =  Y*hm  (19) 

i 

where  x*  is  the  molar  fraction  of  species  i.  Similarly,  the  mixture 
flow  rates  at  inlet  Qin  and  outlet  Qout  can  be  expressed  in  terms  of 
the  ideal  gas  law: 

<*=?£"«  (20) 

i 

After  substituting  Eq.  (17)  and  Eq.  (18)  into  Eq.  (16),  and  integrat¬ 
ing  the  equation,  the  pressure  drop  along  the  channel  was  found  to 
be 

16L 

Pin  ~  Pout  =  —— — 32_[Qin(2Min  +  pout)  +  Qout(/^in  +  ^pout)\  (21) 

The  outflow  pressure  can  be  determined  by  Eq.  (21 ),  and  the  average 
pressure  in  a  segment  can  be  calculated  by 

f^avg  =  2  (Pin  +  ^out)  (22) 

From  discussion  above,  molar  flow  rates  of  every  species  and 
pressure  at  the  segment  outlet  are  determined.  These  properties 
are  used  as  inflow  properties  for  the  subsequent  segment. 


2.2.  Anode/cathode  gas  diffusion  layer  model 


knowing  average  pressure,  we  can  determine  the  partial  pressures 
of  hydrogen  (using  the  anode  GDL  model)  and  oxygen  (cathode  GDL 
model)  as  well  as  water  activity  for  either.  The  hydrogen  and  oxygen 
partial  pressures  are  used  in  the  segment  voltage  model.  Water 
activity  is  used  in  the  membrane  hydration  model  to  determine 
water  transport  through  the  membrane. 

Water  transport  from  the  membrane  to  the  channel  via  GDL  in 
two  forms:  gas  and  liquid;  therefore,  we  need  to  consider  under¬ 
saturated  and  saturated  conditions  separately.  At  under-saturated 
conditions,  water  vapor  transport  direction  depends  on  the  RH  in 
the  channel  and  at  the  membrane/GDL  interface.  At  saturated  con¬ 
dition,  water  generated  in  the  catalyst  layer  will  transport  through 
the  GDL  in  liquid  form.  The  liquid  water  in  the  GDL  not  only  induces 
higher  resistance  to  gas  diffusion  but  also  covers  part  of  the  activa¬ 
tion  sites  on  the  catalyst  layer  and  results  in  reduced  cell  voltage. 

The  Stefan-Maxwell  equation  is  used  to  describe  diffusion 
of  multi-component  gas  mixtures  through  the  GDL  [30].  For  n- 
component  gas  diffusion  through  a  porous  medium,  the  molar 
fraction  gradient  of  species  i,  is  in  the  form: 


VXj  = 


xM  -  X/Nj 

PDfj 


(23) 


where  and  Nj  are  molar  flux  of  species  i  and  j.  PDM  is  the  effec¬ 
tive  pressure  diffusivity  product  of  the  mixture  i-j  in  the  porous 
medium,  and  it  is  related  to  that  in  a  nonporous  system  PDz_j  by 
[31] 

PDf?  =PDi_/e1-5  (24) 

where  £  is  the  porosity  of  the  GDL.  The  pressure  diffusivity  PDt_j 
is  dependent  only  on  temperature  T,  and  can  be  estimated  from 
critical  temperature  Tc,  critical  pressure  Pc  and  molecular  weight 
M  of  components  i  and  j  with  the  following  equation  [30]: 


PDH  =  a 


(PdPcj?/3(TdTcjf/u 


1  1 

Wi  +  Wj 


1/2 


(25) 


Nj  and  Nj  in  Eq.  (23)  can  be  the  molar  flux  of  hydrogen,  oxygen,  or 
water  vapor  through  the  GDL,  and  are  calculated  from  the  segment 
current  and  the  membrane  hydration  model. 


2.2.1.  Under-saturated  condition 

In  the  anode  GDL,  which  contains  hydrogen  and  water  vapor, 
the  water  vapor  molar  fraction  gradient  is  expressed  by  the 
Stefan-Maxwell  equation: 

Jzfdl  =  n - 7T - (^v,gdl%2,gdl  -*H2,gdlNv,gdl)  (26) 

UZ  i  an,avgr-/H2-v 

The  direction  z!  is  defined  in  Fig.  2. 

Since  the  sum  of  the  molar  fractions  of  all  species  is  equal  to  1, 
for  anode  we  have 

*v,gdl  +*H2,gdl  —  1  (27) 

In  Eq.  (26),  Nvgdl  is  the  molar  water  transport  determined  by  the 
membrane  hydration  model,  and  hydrogen  molar  flux  through  the 
GDL  is  equal  to  the  reacted  hydrogen  rate  and  is  calculated  by  Eq. 
(9).  Eq.  (26)  can  be  simplified  by  defining 


Because  GDLs  are  made  of  porous  media,  we  need  to  consider  the 
effect  of  porous  media  on  diffusion  of  gas  mixtures.  Each  species  has 
different  diffusivity,  so  the  molar  fraction  of  species  will  vary  along 
the  diffusion  path.  The  purpose  of  the  GDL  model  is  to  calculate 
the  molar  fraction  at  the  GDL/membrane  interface.  Furthermore,  by 


(28) 

(29) 
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Given  the  boundary  condition  x  =  xv>c/g  at  z!  =  0,  the  water  vapor 
molar  fraction  profile  in  the  GDL  is 

*v,gdi(^  )  —  +  exp(-BiZ  )  ^xv  c/g  —  —  ^  (30) 

The  above  equation  describes  the  distribution  of  water  vapor  molar 
fraction  across  the  anode  GDL  due  to  water  vapor  flux  and  hydrogen 
flux.  At  the  GDL/membrane  interface,  z!  =  tgdi,  the  value  of  the  water 
vapor  molar  fraction  is 

*v,g/m  —  3“  ^XP(— ^lfgdl)  ^v,c/g  —  (31) 

Given  the  water  vapor  molar  fraction  at  the  GDL/membrane 
interface,  we  can  calculate  water  activity  at  the  same  location: 


^an,v,g/m 


*v,g/mBan,avg 
Tv,  sat 


(32) 


Eqs.  (36)  to  (38)  can  then  be  expressed  in  the  matrix  form 


*02,gdl 

*v,gdl 

_ 

-XN2,gdl  _ 

^  xo2,gdi  #3  ~Ba  -B  5 
3_  *v,gdi  =  -B3  B4  -Be 

|_xN2,gdlJ  L°  0  (B5  T  Bq) 

or 

dx 

—  =B  x 

dz 

Eq.  (41)  can  be  solved  by  finding  the  state  transition  matrix 
4>(z)  =  exp(B  •  z) 

and  the  boundary  condition  at  z  =  0: 


(40) 


(41) 


(42) 


x(0)  = 


x02,c/g 

Xv,c/g 
-XN2,c/g  _ 


(43) 


The  hydrogen  partial  pressure  at  the  GDL/membrane  interface 
is  an  important  parameter  to  calculate  segment  voltage  and  it  is 
determined  from 

^an,H2,g/m  —  ^an,avg(l  —  Xv?g/m)  (33) 

If  the  RH  of  gas  flow  in  the  channel  is  different  from  that  at  the 
GDL/channel  interface,  there  will  be  water  vapor  flux  in  between. 
The  molar  flux  of  water  vapor  at  the  GDL/channel  interface  depends 
on  the  inflow  RH  and  is  obtained  from 

Nv,conv  —  Rcon  v(xc/g  —  Xjn  )Aseg  ,conv  (34) 

where  the  convective  mass  transfer  coefficient  Kconv  is  defined  by 
the  Sherwood  number. 

Rconv  =  Sh  cDy  dch  (35) 

where  Dy  is  the  diffusivity  of  species  i  in  the  flow  gas  j.  For  laminar 
flow  and  constant  surface  temperature  conditions  in  a  fuel  cell,  the 
Sherwood  number  Sh  is  constant  and  is  equal  to  3.21  [32]. 

In  the  cathode,  three  species  are  flowing  in  the  channel  and  their 
molar  fraction  gradients  across  the  GDL  are  calculated  from  the 
Stefan-Maxwell  equation: 

4Xo2,gdl  RTseg  fx  02,gdl^v,gdl  —  *v,gdl^02,gdl 

dz  Pc a,avg  y  Do2-v 

+  xQ2,gdl^N2,gdl  ~  *N2,gdlNp2,gdl  \ 

D02-N2  ) 

d*v,gdl  RTseg  /^v,gdl^02,gdl  —  ^02,gdl^v,gdl 

dz  Pc  a,avg  V  °02-v 

^v,gdl^N2,gdl  —  ^N2,gdl^v,gdl 
+  A,-N2 

^XN2,gdl  _  PTseg  /xN2,gdl^v,gdl  —  *v,gdl^N2,gdl 

dz  Pc a,avg  y  T*v-N2 


The  molar  fraction  at  the  GDL/membrane  interface  is  then 
x(tgdi)  =  ^(tgdi)-x(O)  (44) 

Once  the  molar  fractions  of  reactants  at  the  GDL/membrane 
interface  are  determined,  the  partial  pressures  of  reactants  can  be 
calculated  and  used  in  the  segment  voltage  model  to  determine  the 
segment  voltage. 


2.2.2.  Saturated  condition 

When  the  anode  is  saturated,  we  assume  the  vapor  pressures  in 
the  channel  and  in  the  GDL  are  both  equal  to  the  saturated  water 
vapor  pressure.  The  vapor  pressure  is  proportional  to  the  molar 
fraction,  so  the  water  vapor  molar  fraction  in  the  anode  GDL  is  equal 
to  the  saturated  vapor  molar  fraction  and  is  constant: 

*v,gdl(z/)  —  *v,sat  =  5— 1 -  (45) 

1  an,avg 

Since  the  vapor  gradient  is  zero,  the  hydrogen  molar  fraction  is  also 
constant  and  is 


*H2,gdl(z)  —  1  —  Xv,sat 

Hence  the  hydrogen  partial  pressure  is 


(46) 


^H2,g/m  —  *H2,gdl(z/)  ^an,avg  (47) 

Similarly,  when  the  cathode  is  saturated,  water  vapor  molar 
fraction  is  equal  to  the  saturated  molar  fraction  and  is  constant, 
therefore: 


d*v,gdl  =  Q 
dz  - 


(48) 


In  addition  the  sum  of  molar  fractions  of  all  species  is  equal  to 
one 


*02,gdl  +*v,gdl  +*N2,gdl  —  1  (49) 

Given  that  xvgdl  is  constant,  if  we  know  the  nitrogen  molar  fraction, 
the  oxygen  molar  fraction  can  be  determined.  Then  we  use  Eq.  (38) 
to  calculate  the  nitrogen  molar  fraction  from  the  linear  equation: 


|  xN2,gd1^02,gdl  -  Xo2,gdlNN2,gdl  \  pgj 

Do2-N2  ) 

where  the  direction  z  is  defined  in  Fig.  2.  Since  nitrogen  does  not 
react,  NN2  gdl  =  0.  Define 


B3  = 
Be- 


RTseg  gdi  b  _  BTseg  No2,gdl 

Tca,avg  Dq2-v  Tca,avg  Dq2-V 

RTseg  ^v,gdl 

Tca,avg  Dv_n2 


b5- 


RTseg  ^02,gdl 

Pc a,avg  Dq2_n2 


(39) 


4XN2,gdl 

dz 


B5XN2,gdl 


(50) 


Given  the  boundary  condition,  xN2  gdi  =  xN2  c/g  at  z  =  0,  the  nitrogen 
molar  fraction  in  the  cathode  GDL  can  be  expressed  as 


*N2,gdi(z)  =  xN2?c/g  exp(B5z)  (51) 

At  z  =  tgdi,  the  molar  fraction  of  nitrogen  at  the  GDL/membrane 
interface  is  calculated  as 


*N2,g/m  —  xN2,c/g  ^xP(B5fgdl) 


(52) 
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Fig.  3.  Schematic  of  liquid  water  accumulation  in  the  GDL. 


The  partial  pressure  of  oxygen  at  the  same  interface  is  then 


^ca,02,g/m  —  ^ca,avg(l  —  *v,sat  —  *N2,g/m)  (53) 

The  presence  of  liquid  water  in  the  GDL  will  influence  the  dif- 
fusivity  of  gas.  The  effective  diffusivity  Dff  is  dependent  on  the 
porosity  e  and  saturation  s  by 

Df  =  Dif(e)g(s )  (54) 


The  saturation  s  is  the  ratio  of  liquid  water  volume  to  pore  volume 


s  = 


Vw 

Vp 


(55) 


Earlier  studies  [5,10,15,33-35]  suggest  the  influence  of  porosity 
on  diffusivity  to  be  approximated  by  a  polynomial  relationship: 


/(e)  =  e1-5 


(56) 


The  presence  of  liquid  water  reduces  the  diffusion  area  in  the 
GDL  and  its  effect  is  commonly  modeled  by  a  polynomial  function 

g(s)  =  (l-s)m  (57) 


In  this  study  m  =  2  is  used  according  to  Nam  and  Kavinay  [36]. 

It  is  clear  that  liquid  water  accumulation  is  an  important  factor 
that  influences  diffusivities  of  gas  and  cell  performance.  From  neu¬ 
tron  radiography  experiments  [15-17],  water  could  accumulate  in 
the  GDL  under  the  rib.  Thus,  we  propose  a  hypothesis  that  due  to 
the  gas  flow  in  the  channels,  liquid  water  tends  to  move  to  the  GDL 
under  the  rib  and  accumulates  close  to  the  graphite  plate/GDL  inter¬ 
face.  When  the  GDL  cannot  hold  any  more  water,  additional  water 
generated  at  the  membrane  will  push  part  of  the  water  under  the  rib 
into  the  channel,  as  depicted  in  Fig.  3  [28].  In  Ref.  [28],  a  method  to 
differentiate  liquid  water  between  anode  and  cathode  in  the  GDL 
is  developed  based  on  neutron  experimental  data.  From  the  test 
results,  we  observed  a  maximum  water  thickness  tw,max  accumu¬ 
lated  in  the  cathode  GDL  under  the  rib.  In  addition,  large  gas  flow 
rate  in  the  channel  is  also  a  factor  that  influences  water  accumu¬ 
lation  in  the  GDL  under  the  rib.  Thus,  the  average  water  thickness 
(tca.gdi.rib)  in  the  GDL  under  the  rib  in  a  segment  is  approximated  by 


tca,gdl,rib  —  tw,max 


^1  -  exp 


/cdVca!w,gdl,rib  \  \ 
V  (Nca,gas,ch)K  )  ) 


(58) 


where  tw,max  is  the  maximum  water  thickness  that  accumulates 
in  the  GDL  under  the  rib.  This  value  was  determined  by  neutron 
radiography  experiments  [28],  which  is  approximately  50  p,m. 

The  average  liquid  saturation  in  the  GDL  under  the  rib  can  be 
calculated  from 


Sca,gdl,rib 


t, 


ca,gdl,rib  ^seg 


tca,gdl,rib 
tgdl  £ 


(59) 


where  tgdl  is  the  GDL  thickness,  A  is  the  segment  area,  and  e  is  the 
porosity  of  GDL. 

The  liquid  saturation  in  the  GDL  under  the  channel  is  calculated 
from  Refs.  [11,36]: 


KJ  _  PwKKrw  ( dPcap  ^  ds 

c^w,gdl,ch  -  -  Mw/Xw  J  di 


(60) 


At  steady-state,  Nca  wgdl  ch  is  equal  to  the  net  water  flux  from  anode 
to  cathode  and  is  determined  by  the  membrane  hydration  model. 
Then  saturation  s  can  be  obtained  by  solving  Eq.  (60). 

In  this  model,  we  need  to  determine  the  amount  of  liquid 
water  transport  through  the  GDL  under  the  channel  Nca  wgdi  ch  and 
through  the  GDL  under  the  rib  Nca  wgdl  rib.  We  assume  the  water 
transport  through  the  GDL  under  the  channel  is  inversely  propor¬ 
tional  to  the  gas  flow  rate  in  the  channel,  which  is  similar  to  our 
analysis  of  water  accumulation  in  the  GDL  under  the  channel  in 
our  former  study;  i.e.  Cgc=a/Nca  [28], 


^ca,w,gdl,ch  —  at  (61) 

iVca,  react,  ch 

Once  we  determine  the  water  transport  in  the  GDL  under  the 
channel,  we  can  calculate  the  distribution  of  liquid  saturation  in 
the  GDL  by  using  Eq.  (60).  To  simplify  the  problem,  the  capillary 
pressure  gradient  with  respect  to  liquid  saturation  in  Eq.  (60)  is 
assumed  to  be  constant  and  equal  to  22.95  Nm-2  [35].  Then  the 
liquid  saturation  in  the  catalyst  layer  can  be  obtained. 

a,  /3,  and  y  in  Eqs.  (58)  and  (61 )  are  calibrated  based  on  experi¬ 
mental  results.  We  calculate  liquid  saturations  in  the  GDL  under  the 
channel  and  under  the  rib  separately.  The  average  value  between 
these  two  variables  is  then  used  in  the  segment  voltage  model. 


2.3.  Membrane  hydration  model 

The  water  transport  within  membranes  is  represented  by  the 
membrane  hydration  model  shown  in  Fig.  2.  There  are  three  causes 
for  water  flux  in  the  membrane:  electro-osmotic  drag  from  anode 
to  cathode;  back  diffusion  due  to  the  concentration  potential  dif¬ 
ference  between  the  anode  and  cathode;  and  water  generation  at 
the  cathode  catalyst  layer.  These  three  factors  are  explained  in  the 
following. 

The  electro-osmotic  drag  is  defined  as 

^w, osmotic  —  ^osmotic  “p-  (62) 

where  /CosmotjC  is  the  osmotic  drag  coefficient  [3,37-40].  In  this 
study,  we  use  Springer’s  result  [3]: 

Osmotic  =  (63) 

The  water  content  in  the  membrane  Apem,  is  calculated  from  the 
water  activity  of  membrane  apem- 

^pem=0-043-|-17.81CZpem— 39.85dpem-l-36.0apem,  0  <  flpem  <  1 

(64) 


The  average  water  activity  of  anode  and  cathode  is  used  to  cal 
culate  the  water  content  in  the  membrane: 

Qan  +  Qca 


Opem  —  ^ 

In  Eq.  (65),  aan  and  aca  are  the  RH  of  anode  and  cathode. 
The  water  transport  by  back  diffusion  is  expressed  as 


(65) 


i^w,diff  —  ^diff 


ipem 


tgdl  ^seg  £ 


(66) 
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Fig.  4.  Connection  of  six  sub-models  in  the  SIMULINK  environment. 


where  cw  is  the  water  concentration  of  the  membrane  as  defined 
in  Fuller’s  study  [4],  and  tpe m  is  the  membrane  thickness.  Water 
concentration  is  calculated  as 


Cw,an  — 


Ppem 

Mpem 


(67) 


Cw,ca 


Ppem 

Mpem 


A-ca 


(68) 


where  A,an  and  Aca  are  calculated  from  the  water  activity  in  anode 
and  cathode: 

Xan  =  0.043  +  17.81aan-39.85a^n  +  36.0a^n,  0  <  aan  <  1  (69) 

A.ca  =  0.043  +  17.81aca  -  39.85a^a  +  36.0a^a,  0  <  aca  <  1  (70) 


In  Eq.  (66),  Kdiff  is  the  back  diffusion  coefficient  and  is  a  function 
of  temperature  and  water  content  in  the  membrane  [3,41  ]: 

Kji.-Ki  exp  (2416  (jlj  -  -L))  (71) 

where1 


o 

1 

i  o 

?  A-pem  <  2 

10-,O(l+2(Apem-2)) 

,  2  <  A.  pem  <  3 

10-10(3  —  1.167(A.pem  -  3)) 

,  3  <  A,pem  <  4.5 

„  1.25  X  10“10 

?  7 pem  >  4.5 

The  net  water  flux  through  the  anode  GDL  is 

i^an,w,gdl  —  i^w,diff  —  ^w,  osmotic 


(73) 


whereas  that  through  the  cathode  GDL  is 

i^ca,w,gdl  —  osmotic  —  ^w,diff  +  (74) 

where  the  last  term  of  Eq.  ( 74)  is  the  water  generation  at  the  cathode 
catalyst  layer.  The  water  content  is  used  to  calculate  the  membrane 
conductivity  in  the  segment  voltage  model,  and  water  fluxes  in 
anode  and  cathode  are  used  in  the  GDL  models  and  the  flow  channel 
models. 


2.4.  Segment  voltage  model 

The  segment  voltage  model  calculates  voltage  of  each  segment 
at  specific  current  according  to  the  partial  pressures  of  hydrogen 
and  oxygen,  membrane  water  content,  and  temperature.  The  seg¬ 
ment  voltage  can  be  expressed  as 

17seg  =  Wev  —  i7act  —  IA0hm  —  Vconc  (75) 


1  The  value  1.167  is  used  instead  of  1.67.  If  1.67  is  used  in  this  equation,  there  will 
be  a  “discontinuity”  when  calculating  I<x  at  A.pem=4.5  by  using  the  3rd  and  the  4th 
equations. 


where  Vrev,  Vact,  Vohm  and  VConc  are  the  theoretical  reversible  volt¬ 
age,  the  activation  overpotential,  the  ohmic  overpotential,  and  the 
concentration  overpotential,  respectively. 
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The  theoretical  reversible  voltage  is  calculated  from  Ref.  [42]: 
-AG 


Vrev  = 


2  F 

-AC0 
2  F 


RTse g  ln  /  PH2Pg25 


2  F 


h2o 


(76) 


=  1.229  -  0.85  x  10-3(Tseg  -  298.15) 


+  4.3085  x  10-  j  Seg 


ln(^an,H2,g/m)  +  ^  ln(^ca,02,g/m) 


The  partial  pressures  of  hydrogen  Pan,H2,g/m  and  oxygen  Pca,o2,g/m 
come  from  the  anode/cathode  GDL  models. 

The  activation  overpotential  arises  from  the  kinetic  reaction  at 
the  anode  and  cathode.  Due  to  slower  kinetics  of  oxygen  reduction 
at  the  cathode  side,  the  voltage  drop  due  to  activation  overpotential 
is  dominated  by  the  cathode.  The  overpotential  is  modified  from  the 
study  in  Ref.  [25]  as 


Vact  =  0.2  +  0.1[1  -exp(-20iseg)]  (77) 

where  ise g  is  the  current  density  of  a  segment. 

The  ohmic  overpotential  is  due  to  the  internal  resistance  of  a 
segment  and  is  expressed  as 


^ohm  —  fsegPseg  (78) 

where  /seg  is  the  segment  current.  The  resistance  of  the  segment  is 
the  sum  of  all  components  through  which  current  flows  and  contact 
resistance.  These  components  are  membrane,  GDLs,  and  flow  field 
plates: 


Pseg  —  Ppem  +  Pgdl  +  Opiate  +  ^contact  (79) 

The  conductivities  of  GDL  and  graphite  plates  are  typically  much 
larger  than  that  of  the  membrane,  so  it  is  not  necessary  to  consider 
their  resistances.  Thus  only  the  membrane  resistance  and  contact 
resistance  are  considered  in  this  model.  The  membrane  resistance 
is  obtained  from 


D  _  fpem 

rtpem  —  - r 

O+em^seg 


(80) 


where  the  membrane  conductivity  crpem  is  a  function  of  tempera¬ 
ture  and  water  content  in  the  membrane,  and  is  expressed  in  the 
form  [3]: 


O+em  —  (bnApem  —  b\2)exp 


bi3 


(81) 


where  b\\ ,  b\2,  and  bu  are  empirically  determined  from  our  exper¬ 
imental  results. 

Concentration  overpotential  results  from  the  change  in  concen¬ 
tration  of  the  reactants  as  they  are  consumed  in  the  reaction.  The 
concentration  overpotential  derived  from  the  Nernst  equation  is 
modified  from  Ref.  [43].  In  addition,  the  flooding  effect  should  also 
be  considered,  which  reduces  the  activation  area  of  the  catalyst, 
so  the  maximum  current  density  is  reduced  when  liquid  water 
appears  in  the  catalyst  layer.  The  modified  concentration  overpo¬ 
tential  is  expressed  as 


Vconc  -  b2 1 42i  In  ( -  iuJjf-s))  (82) 

where  b2\  and  b2 2  are  coefficients  to  be  determined  by  experimen¬ 
tal  data. 

In  the  above  calculations,  the  current  density  in  each  segment 
is  assumed  to  be  known  and  the  same.  However,  the  cell  voltage  of 
all  segments  should  be  the  same,  and  the  difference  in  humidity, 
reactant  pressure,  etc.  resulted  in  different  current  density.  After 
each  segment  voltage  model  calculates  its  voltage,  actual  segment 


c 


Yes 

♦ 


End  Vcctt=  Vm 


j 


Fig.  5.  Flow  diagram  of  solving  cell  voltage. 


current  can  be  corrected  by  enforcing  all  the  segment  voltages  to 
be  the  same. 

2.5.  Cell  voltage  calculation 

The  models  presented  in  the  previous  sub-sections  are  imple¬ 
mented  in  the  SIMULINK  environment.  The  block  diagrams  of  six 
sub-models  of  a  segment  are  shown  in  Fig.  4.  Given  the  inflow  prop¬ 
erties  and  (initial  guess  of)  segment  current,  the  segment  model 
calculates  the  segment  voltage.  Since  the  current  density  is  actu¬ 
ally  not  uniform  throughout  the  active  area  of  the  cell.  The  cell 
voltage  is  determined  iteratively  by  the  process  shown  in  Fig.  5.  At 
the  beginning,  we  guess  an  inlet  gas  pressure  Pin  and  a  cell  current 
/ceii.  The  initial  guess  for  segment  currents  /seg  2-,  is  set  to  be  one  15th 


Table  1 

Specification  of  the  single  cell  used  in  this  study 


Parameter 

Value 

Cell  active  area  (Acen) 

100  cm2 

Channel  depth  (Hch) 

1  mm 

Channel  width  (Wch) 

1.6  mm 

Rib  width  (Wrib) 

1.7  mm 

Anode  channel  number  in  a  segment  (Zan) 

6 

Cathode  channel  number  in  a  segment  (Zca) 

10 

Anode  channel  length  in  a  segment  (Lan) 

3.33  cm 

Cathode  channel  length  in  a  segment  (Lca) 

2  cm 

GDL  thickness  (tgd  1) 

184  (Jim 

GDL  porosity  (e) 

0.725 

Dry  membrane  thickness  (tpem) 

25  (Jim 

Dry  membrane  density  (ppem) 

2000  kg  irr3 

Dry  membrane  equivalent  weight  (Mpem) 

1.1  kg  mol-1 

Table  2 

Parameter  values  that  obtained  from  literature 
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Parameter 

Value 

Absolute  permeability  of  GDL 

00 

1  x  10“8  m-2  [11] 

Relative  permeability  of  GDL 

(Krw) 

S  [12] 

Dynamic  viscosity  of  hydrogen 
(Mh2) 

9.5493  x  10-6  Nsm-2  [29] 

Dynamic  viscosity  of  oxygen 

(Mo2) 

2.2379  xlO-5Nsm-2  [29] 

Dynamic  viscosity  of  nitrogen 

(Mn2) 

1.9260  x  10“5  Nsm”2  [29] 

Dynamic  viscosity  of  water 
vapor  (/xv) 

4.6657  xl0“4  Nsm-2  [29] 

Pressure-diffusivity  product  of 
water  vapor  and  hydrogen 
pair  (PDh2-v) 

16.6801  Pa  m2  s_1  [30] 

Pressure-diffusivity  product  of 
water  vapor  and  oxygen  pair 

(PDo2-v) 

3.2890  Pa  m2  s”1  [30] 

Pressure-diffusivity  product  of 
water  vapor  and  nitrogen 
pair  (PDv_n2  ) 

3.4400  Pa  m2s~1  [30] 

Pressure-diffusivity  product  of 
oxygen  and  nitrogen  pair 
(PDo2_n2  ) 

2.5504  Pa  m2s-1  [30] 

Diffusivity  of  water  vapor  in 
hydrogen  (Dv_h2  ) 

9.3940  x  10“3  m2  s”1  [30] 

Diffusivity  of  water  vapor  in  air 

(Dv- air) 

2.6560  x  10-3  m2  s”1  [30] 

of  the  cell  current,  Ice\\.  The  segment  voltage  Vseg  i  is  determined  by 
the  segment  model  described  in  the  previous  section. 

If  the  difference  between  the  maximum  and  minimum  segment 
voltages  is  not  within  an  acceptable  range,  the  segment  currents 
need  to  be  corrected.  Based  on  the  typical  polarization  curve  of  a 
fuel  cell,  the  segment  with  higher  voltage  should  increase  its  cur¬ 
rent,  and  that  with  lower  voltage  should  decrease  its  current.  To 
increase  the  iteration  speed,  the  increase  in  segment  current  is  set 
to  be  proportional  to  the  voltage  difference  while  keeping  the  cell 
current  constant.  In  addition,  the  inlet  pressures  of  the  anode  and 
cathode  are  also  adjusted  to  keep  the  outlet  pressure  the  same  as 
the  ambient  pressure. 

2.6.  Tuning  parameters 

We  designed  a  single  cell  and  conducted  neutron  radiography 
experiments  to  adjust  water  related  parameters  in  this  model.  The 
fuel  cell  specification  and  corresponding  parameters  are  listed  in 
Table  1 .  Those  parameter  values  obtained  from  the  literature  are 
listed  in  Table  2.  The  parameters  that  are  adjusted  in  this  paper  are 
listed  in  Table  3. 

In  the  ohmic  overpotential  region,  membrane  conductivity  is 
mainly  affected  by  water  content  in  the  membrane,  as  shown 


Fig.  6.  Comparison  of  experimental  results  and  modeling  results  for  different  cath¬ 
ode  inlet  RH. 

in  Eq.  (81).  Thus,  bn ,  b12,  and  b\2  in  Eq.  (81)  and  the  contact 
resistance  ^contact  in  Eq.  (79)  were  tuned  to  match  the  calculated 
I-V  curve  to  with  experimental  data  in  the  ohmic  overpoten¬ 
tial  region,  a  and  y  in  Eq.  (58)  and  f3  in  Eq.  (61)  were  adjusted 
based  on  water  thickness  obtained  by  neutron  radiography  exper¬ 
iments.  Once  a ,  /3,  and  y  were  determined,  liquid  saturation  can 
be  determined.  Subsequently,  b2 1,  b22 ,  and  ilimit  in  Eq.  (82)  are 
adjusted  according  to  the  concentration  overpotential  region  of 
I-V  curve. 

3.  Results  and  discussions 

In  practical  applications,  due  to  space  and  cost  considera¬ 
tions,  it  is  common  to  humidify  only  the  cathode  reactant.  Thus, 
we  will  focus  on  the  influence  of  RH  of  cathode  inlet  on  cell 
performance  and  water  transport.  Fig.  6  compares  the  cell  per¬ 
formance  curves  obtained  by  this  model  and  by  experiments.  The 
cell  operating  at  low  cathode  inlet  RH  of  50%  shows  lower  cell 
performance.  This  is  because  the  under-saturated  air  takes  water 
from  the  membrane,  resulting  in  low  membrane  hydration  and 
conductivity. 

In  this  study,  each  segment  is  regarded  as  a  lumped  model. 
Thus,  at  the  end  of  iteration  procedure,  only  average  values  will 
be  obtained.  For  visual  aids,  colorful  pictures  are  created  through 
interpolation  and  extrapolation  of  these  values,  which  are  marked 
at  the  center  of  each  segment. 

3.1.  Distribution  of  current  density  and  water  content  in  the 
membrane 


Table  3 

Parameters  that  tuned  based  on  experimental  data 


Parameter 

Eq.  number 

Tuned  value 

^contact 

(79) 

0.047 

a 

(58) 

-0.012 

P 

(61) 

1  x  10-9 

Y 

(58) 

2 

bn 

(81) 

0.195 

bn 

(81) 

0.326 

bn 

(81) 

350 

£>21 

(82) 

-0.75 

b22 

(82) 

7 

Mimit 

(82) 

1.2 

Membrane  dehydration  increases  ohmic  overpotential,  and 
could  even  cause  irreversible  damage  to  the  membrane.  Since 
the  membrane  conductivity  dominates  (conductivity  of  mem¬ 
brane:  ~2.75  m_1 ;  GDL:  ~1250  m-1 ;  graphite  plate: 

Table  4 

List  of  selected  operating  conditions 
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Current  density 

distribution 

(%) 


I 


(a) 


Fig.  7.  Distribution  of  (a)  and  (b)  current  density;  (c)  and  (d)  water  content  in  the  membrane;  (e)  and  (f)  RH  in  the  anode  channel;  (g)  and  (h)  RH  in  the  cathode  channel.  Left 
figures:  case  1.  Right  figures:  case  2. 


~1  x  105  m_1),  current  density  distribution  is  highly  depen¬ 

dent  on  the  water  content  in  the  membrane,  which  in  terms 
significantly  influences  the  cell  performance  and  reliability. 

Figs.  7(a)  and  (b)  and  8(a)  and  (b)  show  current  density  distri¬ 
bution  of  four  selected  operating  conditions,  as  listed  in  Table  4.  All 
figures  suggest  that  maximum  current  are  located  near  the  anode 


outlet,  the  region  with  maximum  water  content.  By  comparing 
Fig.  7(a)  and  (b)  (high  current  density  cases),  we  observe  that  at 
fully  humidified  cathode  condition,  the  maximum  local  current  is 
approximately  1.4  times  the  minimum  local  current.  However,  this 
ratio  increases  to  2  at  low  cathode  humidity  condition.  In  Fig.  8(a) 
and  (b),  for  low  current  density  cases,  the  maximum  local  current 


Y.-S.  Chen,  H.  Peng  /  Journal  of  Power  Sources  185  (2008)  1179-1192 


1189 


Current  density 

distribution 

(%) 


Case  3 


(a) 


Case  4 


M 


2 

0 


Water  content 


Anode  RH 

(%) 


(C) 


Cathode  RH 

(%) 


(e) 


(g) 


B 

I 


+  100.00  +  100.00 

+  100.00  +  100.00 

+  100.00  +  100.00 

+  100.00  +  100.00 

+  100.00  +  100.00 


A 


■  100. 

-  worn 

•  ico  bo 


+  100 

+  100 


•  — — 


+  50.12 

+  50.11 

+  50.1  M 

+  50.37 

+  50.36 

+  50.6  > 

+  52.50 

+  52.26 

+  52  2  1 

+  56.21 

+  56.15 

+  55 

+  60  40 

+  60.57 

Fig.  8.  Distribution  of  (a)  and  (b)  current  density;  (c)  and  (d)  water  content  in  the  membrane;  (e)  and  (f)  RH  in  the  anode  channel;  (g)  and  (h)  RH  in  the  cathode  channel.  Left 
figures:  case  3.  Right  figures:  case  4. 


is  1.5  and  1.6  times  of  the  minimum,  respectively.  The  result  sug¬ 
gests  that  at  low  current  density,  cathode  inlet  RH  has  less  influence 
on  current  density  distribution.  Low  cathode  humidity  causes  low 
water  content,  as  shown  in  Fig.  8(d). 

A  few  methods  to  measure  current  density  distribution  were 
developed  in  the  literature  [44-47].  For  example,  divided  current 


collectors  can  be  used  to  measure  local  current  flow.  Yoshioka  et 
al.  [44]  compared  the  distribution  of  current  densities  at  different 
inlet  gas  RH  levels.  Their  results  show  that  the  region  close  to  air 
inlet  has  lower  current  density.  The  trend  is  more  significant  with 
dry  inlet  air.  Liu  et  al.  [45]  measured  current  density  distribution 
of  a  fuel  cell  with  one  serpentine  flow  channel.  Their  findings  are 
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Fig.  9.  Quantity  of  water  transport  (mol  s-1 )  across  the  MEA.  (a)  Case  1 ;  (b)  case  2;  (c)  case  3;  (d)  case  4. 
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the  same  as  [44]  qualitatively,  except  that  with  fully  humidified  air, 
current  density  of  the  area  near  outlet  is  lower  due  to  flooding.  In 
our  study,  the  effect  of  flooding  on  reducing  current  density  near 
outlet  is  not  obvious.  A  possible  reason  is  that  anode  reactant  is  not 
humidified  so  the  excess  water  transports  to  the  anode.  Another 
reason  is  that  we  use  straight  parallel  cathode  channels  so  liquid 
water  is  quickly  removed  from  channels. 


3.2.  Distribution  ofRH  in  the  flow  channel  and  water  transport  in 
the  MEA 

Figs.  7(e)  and  (f)  and  8(e)  and  (f)  show  the  RH  distribution  in 
the  anode  flow  channel.  For  all  cases,  RFI  increases  along  the  flow 
direction  because  dry  hydrogen  gradually  uptakes  water  vapor  that 
comes  from  the  cathode  by  back  diffusion.  Fig.  9  compares  water 
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Fig.  10.  Comparison  of  modeling  results  and  experimental  data,  (a)  Distribution  of  liquid  saturation  in  the  GDL  under  the  rib  when  cathode  inlet  RH  =  100%;  (b)  distribution 
of  liquid  saturation  in  the  GDL  under  the  channel  when  cathode  inlet  RH  =  100%;  (c)  distribution  of  liquid  saturation  in  the  GDL  under  the  rib  when  cathode  inlet  RH  =  50%; 
(d)  distribution  of  liquid  saturation  in  the  GDL  under  the  channel  when  cathode  inlet  RH  =  50%. 


transport  through  the  membrane  and  through  both  anode  and  cath¬ 
ode  GDLs.  The  magnitude  indicates  the  quantity  of  water  transport 
and  positive  represents  the  direction  from  the  anode  to  the  cathode. 
At  high  current  density,  water  transport  from  the  anode  channel  to 
the  membrane  through  the  anode  GDL  was  observed  in  the  down¬ 
stream  segments,  as  shown  in  Fig.  9(a)  and  (b).  The  high  current 
density  in  those  segments  results  in  high  electro-osmotic  drag  from 
the  anode  to  the  cathode,  which  is  stronger  than  the  back  diffusion 
from  the  cathode  to  the  anode;  hence,  water  in  the  anode  reac¬ 
tant  supplies  for  the  difference.  However,  as  shown  in  Fig.  7(a),  the 
anode  RH  does  not  decrease  along  flow  channel.  Because  of  hydro¬ 
gen  consumption  along  flow  channel,  the  molar  fraction  of  water 
vapor  increases  along  flow  channel.  That  also  explains  why  RH  in 
the  anode  downstream  segments  barely  increased. 

Fig.  9(c)  and  (d)  shows  that  the  electro-osmotic  drag  of  the  seg¬ 
ments  near  outlet  in  case  3  is  significantly  higher  than  that  in  case  4, 
although  current  densities  of  both  cases  are  quite  similar,  as  shown 
in  Fig.  8(a)  and  (b).  According  to  Eq.  (63),  the  electro-osmotic  drag 
coefficient  is  a  function  of  water  content.  Thus,  high  water  content 
in  the  membrane  is  attributed  to  humidified  cathode  reactant. 

Figs.  7(g)  and  8(g)  show  fully  saturated  reactant  in  the  cathode 
channel  throughout  the  active  area.  Since  cathode  inlet  reactant 
is  fully  humidified,  the  generated  water  transports  to  the  chan¬ 


nel  in  liquid  form.  In  the  cases  of  under-saturated  cathode  inlet 
reactant,  the  cathode  RH  gradually  increases  along  the  gas  flow 
direction,  as  shown  in  Figs.  7(h)  and  8(h).  The  increase  is  due  to 
water  generation  and  also  oxygen  consumption  along  the  flow 
channel,  resulting  in  increased  molar  fraction  of  water  vapor.  Water 
transports  through  GDLs  influences  diffusivities  of  reactant  and 
may  also  form  liquid  water  in  the  GDL,  which  increases  concen¬ 
tration  overpotential. 

3.3.  Distribution  of  water  accumulation 

Cathode  liquid  saturation  in  the  GDL  affects  cell  performance 
because  liquid  water  may  cover  the  reaction  sites  in  the  catalyst 
layer  or  block  the  pathway  of  gas  flow  through  the  GDL.  The  liq¬ 
uid  saturation  in  the  cathode  GDL  is  affected  by  many  factors:  cell 
temperature,  permeability  and  hydrophobicity  of  GDL,  net  water 
flux  through  the  GDL,  and  cathode  flow  RH.  In  this  model,  temper¬ 
ature,  permeability  and  hydrophobicity  of  GDL  are  assumed  to  be 
constant.  Under-saturated  condition,  net  water  flux  through  GDL 
is  determined  by  operating  current  density.  Higher  current  density 
results  in  higher  electro-osmotic  drag  and  water  generation,  which 
causes  higher  net  water  flux  through  GDL,  as  shown  in  Fig.  9(a)  and 
(b). 
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Fig.  10  compares  average  liquid  saturation  in  the  cathode  GDL 
of  modeling  results  and  experimental  data.  The  15  subplots  corre¬ 
spond  to  15  segments  of  the  active  area.  Fig.  10(a)  and  (c)  suggests 
that  liquid  water  accumulate  in  the  GDL  under  the  rib.  Results 
from  other  neutron  radiography  experiments  also  show  this  phe¬ 
nomenon  [48,49].  Slight  liquid  saturation  in  the  GDL  under  the  rib 
is  observed  at  the  segment  near  cathode  inlet  and  increases  toward 
the  outlet.  Maximum  liquid  saturation  is  approximately  0.4.  The 
low  water  accumulation  in  the  cathode  GDL  at  high  current  den¬ 
sity  is  likely  due  to  the  high  gas  flow  rate  at  high  current  density. 
Fig.  10(c)  shows  some  liquid  saturation  in  the  GDL  under  the  rib 
when  cathode  inlet  RH  is  50%.  However,  there  is  almost  no  liquid 
water  in  the  GDL  under  the  channel,  as  shown  in  Fig.  10(d).  This  is 
because  the  gas  flow  in  the  flow  channel  is  under-saturated. 

4.  Conclusions 

In  this  study  a  segmented  fuel  cell  model  was  developed.  Each 
segment  is  viewed  as  a  lumped  model,  which  consists  of  six  inter¬ 
acting  sub-models,  and  is  connected  based  on  the  reactant  flow 
directions.  This  model  is  calibrated  based  on  experimental  results 
and  is  able  to  describe  liquid  water  saturation  in  the  GDL  under  the 
channel  and  in  the  GDL  under  the  rib.  The  calibrated  model  was 
used  to  investigate  distributions  of  current  density,  water  content 
in  the  membrane,  RH  in  the  anode/cathode,  and  water  accumula¬ 
tion  in  the  GDL  under  the  channel/rib  as  well  as  water  transport  in 
the  MEA. 

Modeling  results  show  that  cathode  inlet  RH  has  significant 
influence  on  the  uniformity  of  water  content  in  the  membrane  and 
current  density.  At  low  cathode  inlet  humidity  and  high  load,  the 
maximum  local  current  density  is  twice  that  of  the  minimum  local 
current  density.  Cathode  humidity  has  less  influence  on  the  uni¬ 
formity  of  current  density  at  low  current  load.  In  this  study,  the 
influence  of  water  accumulation  on  current  density  is  not  obvious. 
That  could  be  due  to  un-humidified  hydrogen,  short  flow  channels, 
and  operating  current  density. 

Water  transport  mechanisms  across  the  MEA  were  demon¬ 
strated  in  this  study.  The  amount  of  water  transport  in  different 
location  was  predicted  by  this  segmented  model.  The  results  show 
current  density  and  the  amount  of  water  transport  influence  with 
each  other.  This  result  also  provides  useful  information  in  plac¬ 
ing  the  inlet/oulet  of  anode/cathode  when  we  design  flow  field 
patterns. 

A  model  describing  water  accumulation  in  the  GDL  was  pro¬ 
posed  in  this  study.  Liquid  water  tends  to  accumulate  in  the  GDL 
under  the  rib  due  to  the  suppression  of  gas  flow  in  the  channel. 
Maximum  liquid  saturation  in  the  GDL  under  the  rib  is  observed 
near  the  outlet  and  the  value  is  approximately  0.4  in  this  flow  field 
design. 

This  study  investigated  the  water  content  in  the  membrane  and 
liquid  water  accumulation  in  the  GDL  for  a  specific  designed  single 
cell.  For  the  future  study,  the  concept  of  modeling  will  be  applied 
on  different  flow  field  designs  to  study  the  influence  of  flow  field 
design  on  water  distribution  in  a  fuel  cell. 
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